LISA detections of massive black hole inspirals: parameter extraction errors due to 

inaccurate template waveforms 
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The planned Laser Interferometer Space Antenna (LISA) is expected to detect the inspiral and 
merger of massive black hole binaries (MBHBs) at 2 < 5 with signal-to-noise ratios (SNRs) of 
hundreds to thousands. Because of these high SNRs, and because these SNRs accrete over periods 
of weeks to months, it should be possible to extract the physical parameters of these systems with 
high accuracy; for instance, for a ~ 1O^M0 MBHBs at z = 1 it should be possible to determine 
the two masses to ~ 0.1% and the sky location to ~ 1°. However, those are just the errors due to 
noise: there will be additional "theoretical" errors due to inaccuracies in our best model waveforms, 
which are still only approximate. The goal of this paper is to estimate the typical magnitude of 
these theoretical errors. We develop mathematical tools for this purpose, and apply them to a 
somewhat simplified version of the MBHB problem, in which we consider just the inspiral part of 
the waveform and neglect spin-induced precession, eccentricity, and PN amplitude corrections. For 
this simplified version, we estimate that theoretical uncertainties in sky position will typically be 
~ 1°, i.e., comparable to the statistical uncertainty. For the mass and spin parameters, our results 
suggest that while theoretical errors will be rather small absolutely, they could still dominate over 
statistical errors (by roughly an order of magnitude) for the strongest sources. The tools developed 
here should be useful for estimating the magnitude of theoretical errors in many other problems in 
gravitational- wave astronomy. 

PACS numbers: 04.25.Nx,04.30.Db,04.80.Nn,95.75.Wx,95.85.Sz 



I. INTRODUCTION 

The inspirals and mergers of massive (~ \Q^Mq) 
black-hole binaries (MBHBs) are likely to be the 
strongest gravitational-wave (GW) sources detected by 
LISA, the planned Laser Interferometer Space Antenna. 
MBHBs at redshifts 2 < 5 will be detected by LISA with 
optimal matched-filtering signal-to-noise ratios (SNRs) 
~ 10^-10'^. Because of these high SNRs, and because 
these sources emit GWs in the LISA band for weeks to 
months, it should be possible to determine the source pa- 
rameters to very high accuracy. For instance, from the 
inspiral waveform of MBHB systems at z = 1 it should 
be possible to extract the two masses to within ~ 0.1%, 
the spins to within ~ 0.1-1%, the luminosity distance Dl 
to ~ 0.1-1%, and the sky location to within - 0.1-1° [l|. 

More precisely, these are the sizes of the statistical er- 
rors due to random noise, as calculated using the Fisher- 
matrix formalism and (over the years) increasingly faith- 
ful models {templates) oi the gravitational waveforms 
(see, e.g., Refs. [U, H, 1, i; see also Ref. @ about the 
conditions required for the Fisher-matrix formalism to be 
applicable). While waveform models have improved, they 
are still only approximate versions of the true general- 
relativistic (GR) waveforms, and using them to fit the 
detector data will incur an additional theoretical error in 
estimating source parameters. Note that, to lowest order, 
while the statistical error due to noise scales roughly as 
1/SNR, the theoretical error will be independent of SNR. 
For precisely this reason, we might well worry that for the 
highest SNR sources - and especially MBHB inspirals 
- theoretical errors could dominate the total parameter- 



estimation error. This basic point has already been made 
by Berti [g]; however we believe that this paper provides 
the first fairly realistic estimates of the likely magnitude 
of such theoretical errors. 

In recent years numerical relativity (NR) has made 
great strides toward solving Einstein's equations for the 
final stages of MBHB inspirals, including the merger and 
ringdown phases [7, 's', and work has begun on the 
comparison of NR waveforms with the last cycles of post- 
Newtonian (PN) waveforms [lOjilUi- Still, the PN expan- 
sion is still the best available method for calculating the 
vast majority of the ~ lO^'^-lO^ GW cycles observed by 
LISA, and it is likely to remain so for many years. In 
the PN approximation, the GR solution is written as a 
"Newtonian" solution plus corrections given as a power 
series in d/c [l^ . Our convention is to refer to correc- 
tions of order {v/cY to either the instantaneous orbital 
elements or their decay rates as fcPN corrections. Thus 
we refer to the lowest-order waveform, corresponding to 
a quasi-Newtonian orbit that is slowly decaying accord- 
ing to the energy-loss rate given by the quadrupole for- 
mula, as the OPN or Newtonian waveform. For circu- 
lar nonspinning BHs, approx_imate waveforms have been 
calculated to 3.5PN order [lj|. For BHs with signifi- 
cant spins (which cause the orbital plane to precess if 
the spin axes are not exactly aligned with the binary's 
orbital angular momentum) , the state of the art is some- 
what less advanced, as the corresponding waveforms have 
been calculated through 2.5PN order [lJ,[lB|- Unfortu- 
nately there is as yet no general algorithm for extending 
PN calculations to arbitrarily high order, and each new 
order appears to require much additional effort. 

We were particularly motivated in researching this pa- 
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per by the question of how much LISA's angular resolu- 
tion is limited by the accuracy of the best available PN 
templates. This question is important because feasible 
searches for electromagnetic counterparts to gravitation- 
ally observed MBHB mergers are likely to require GW- 
derived error boxes no larger than a few square degrees 
on the sky. The statistical errors in sky position will 
generally satisfy this, but what about theoretical errors? 
Put another way: given the great effort and cost of build- 
ing LISA, it seems that a reasonable goal for waveform- 
modeling theorists is that parameter-estimation errors 
due to inaccurate templates should be no larger than 
those arising from noise. This paper develops tools that 
will help theorists estimate how close they currently are 
to this goal, where future effort might best be placed, 
and when they can reasonably stop. 

At this point, the reader might well be wondering: if 
full waveforms /igr from NR are not yet available, by 
what comparison can we judge the accuracy of high-order 
PN solutions? Strictly speaking, of course we cannot; 
nevertheless, there at least two ways that we can explore 
this problem. First, in the limit of very small mass ra- 
tios {M2/M1 <^ 1), the smaller BH can be treated as 
a perturbation on the spacetime of the larger BH, and 
the problem can be solved up to order {M2IM1Y using 
BH perturbation theory. This is straightforward when 
both BHs are nonspinning, so at least in this case we 
can compare the high-order PN solution to a solution of 
(practically) arbitrary accuracy. 

Second, since the PN expansion is known to converge 
rather slowly for the late stages of inspiral [H, [iTj . 
we might estimate that the waveform error ft.GR(i) ~ 
/i3.5PN(i) is comparable in size to /i3.5PN(i) — /i3.0PN(i). 
The latter observation leads to a particularly simple 
method for estimating the magnitude of theoretical er- 
rors, which we implement in this paper. In effect, we 
"replace" the true waveform /iGR(i) by /i3.5PN(i) and its 
best approximation /i3.5PN(i) by h^.ovN{t), and directly 
calculate the resulting errors. Indeed, essentially this 
same starting point was taken by Canitrot p^ . who in- 
vestigated theoretical errors in the context of stellar-mass 
compact binaries detected by ground-based GW detec- 
tors. Canitrot found, e.g., that fitting 2.0PN templates 
to 2.5PN waveforms (using the frequency-weighting ap- 
propriate for the Virgo detector noise curve) leads to the- 
oretical errors in the chirp mass of order ~ 1% (though 
small, this is larger than the estimated statistical error), 
and errors in the symmetric mass ratio of order ~ 50%. 

While we believe that this method leads to a useful first 
estimate of theoretical errors, it ignores data-analysis 
strategies that might be used to mitigate them. In this 
paper, we investigate two such methods of reducing er- 
ror. The first involves deweighting the waveform in the 
last stages of inspiral, where the PN approximation (at 
any given order) is most uncertain. The second strat- 
egy involves using as templates hybrid waveforms that 
include all available test-mass PN corrections, which are 
known to much higher PN order than the corresponding 



comparable-mass terms. 

The plan of this paper is as follows. In Sec.|TT]we give 
a rather general geometrical description of parameter- 
estimation errors in GW data analysis, and we develop 
two tools for estimating the sizes of theoretical errors. 
The first tool relies on the integration of an ordinary dif- 
ferential equation (ODE) to find the approximate tem- 
plate that best fits a given true GW signal (and vice 
versa). The second tool is a one-step substitute for the 
ODE, which is generally accurate enough for our pur- 
poses and considerably faster to evaluate. We feel that 
both of these tools could be profitably applied to many 
similar problems, especially in ground-based GW astron- 
omy. In Sec. mil we describe the PN waveforms that we 
adopt as our model waveforms for MBHB inspirals. Since 
in this paper we attempt only a first pass at the full prob- 
lem, for simplicity a) we employ the restricted PN ap- 
proximation, which retains high-order PN corrections to 
the phase of the waveform's quadrupole component, but 
neglects both higher-multipole components and PN cor- 
rections to the waveform's amplitude; b) we neglect the 
spin-induced precession of the orbital plane; and c) we 
neglect the (post-inspiral) merger and ringdown phases 
of the waveform. It should be straightforward (if labo- 
rious) to restore these neglected pieces in future work. 
In Sec, mil we also summarize our treatment of the LISA 
response function, which relies on a low-frequency ap- 
proximation 2] that is quite adequate for our present 
purposes. Our results are presented in Sec. IIVI In Sec. 
IVl we summarize our conclusions and list some directions 
for future research. 

Throughout this paper, we use geometrical units in 
which G = c ~ I. Therefore everything can be measured 
in a fundamental unit of seconds. For familiarity, we 
sometimes express quantities in terms of yr, Mpc, or M©, 
which are related to our fundamental unit by I yr = 
3.1556 X 10^ s, 1 Mpc = 1.029 x s, and IM© = 
4.926 X 10-6 s. 



II. GENERAL FORMALISM FOR ESTIMATING 
STATISTICAL AND SYSTEMATIC ERRORS 

A. A review of signal analysis 

This section briefly reviews the basic formulas of signal 
analysis, partly to fix notation. For a more complete 
discussion, see Refs. [l9| or [20j . 

The output of D detectors can be represented by the 
vector s"(t), a = 1, 2, . . . , I?. It is often convenient to 
work with the Fourier transform of the signal; our con- 
vention is 

/oc 
e^^'^'s"{t)dt. (1) 
-00 

The detector output s" (t) is the sum of GWs ft." (t) plus 
instrument noise (t) . We assume that the noise is sta- 
tionary and Gaussian; here "stationarity" means that the 
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different Fourier components n°'{f) of the noise are un- 
correlated; thus we have 

{n^{f)n^{n)^\Kf-nsf{f). (2) 

where "( )" denotes averaging over reahzations of the 
noise, and S'j^^{f) is the (one-sided) spectral density of 
the noise. By contrast, "Gaussianity" imphes that each 
Fourier component has Gaussian probabihty distribu- 
tion. Under these assumptions, we obtain a natural inner 
product on the vector space of signals: given two signals 
g^it) and k°'{t), we define (g | k) by 

(g|k) = 4Re/ [Su[f)-'Ur{fr~kP{f)df, (3) 
Jo 

where we sum over repeated indices. Here and below 
we reduce index clutter by using boldface to represent 
multiple-detector abstract signal vectors (e.g., k will de- 
note both fc"(t) and k°'{f), which are the same vector 
represented in different bases). 

In terms of the inner product ([3]), the probability for 
the noise to have some realization no is just 



p{n = no) cx 



.-(no |no)/2 



(4) 



Thus, if the actual incident waveform is h, the probability 
of measuring a signal s in the detector output is propor- 
tional to exp{— (s— h I s— h)/2}. Correspondingly, given a 
measured signal s, the gravitational waveform h that best 
fits the data (in the maximum-likelihood sense) is the one 
that minimizes the quantity (s — h | s — h) . This crite- 
rion has a simple geometric interpretation. If N is the 
number of source parameters, the parametrized space of 
waveforms {h(0*)} is an TV-dimensional manifold embed- 
ded in the vector space of all possible measured signals. 
Given the measured signal s, the best-fit waveform h(6'^f ) 
is the point on the waveform manifold that lies closest to 
s according to the distance (s — h(0^f) | s — h(6'^f)). The 
vector from h(6'f,f) to s is then normal to the waveform 
manifold at dl^f. It is also easily shown that the (matched- 
filtering) SNR of h is just its norm (h|h)^/^. 

For a given incident GW, different realizations of the 
noise will give rise to somewhat different best-fit parame- 
ters. For large SNR, however, the best-fit parameters will 
assume a normal distribution centered on the correct val- 
ues. Specifically, if 91^ are the true physical parameters, 
and 0bf(n) = 91^ + A9^{n) are the best-fit parameters. 



then for large SNR the parameter-estimation errors A9^ 
have the normal distribution 

p(A0') = AAe-^^-^™. (5) 

Here Fy is the Fisher information matrix defined by 



dh 



dh 



(6) 



and J\f — ^ det(r/27r) is the appropriate normalization 
factor. For large SNR, the variance-covariance matrix of 
the errors A6'* is given by 

{M'M') = {T'^y^ + O(SNR-i) (7) 

(see, e.g., Ref. [1]). 

B. Theoretical parameter-estimation errors 

We now give a simple geometric picture of theoreti- 
cal errors due to inaccurate templates. Consider the fol- 
lowing two manifolds embedded in the vector space of 
possible data streams: a manifold of the true GR wave- 
forms {h.QY{{9^Y\ , and a different manifold of approximate 
waveforms {hAp(^*)} generated by some approximation 
procedure. Since the two manifolds can be parametrized 
by the same physical coordinates (the masses of the bod- 
ies, their spin vectors, etc.), there is a natural one-to-one 
mapping between waveforms of the same coordinates on 
the two manifolds. 

Now consider the measured detector data s = 
hGR(^tr) + where n is again the detector noise. Since 
in practice we have access only to the approximate wave- 
forms, the hAp(0bf) that lies closest to the data s yields 
our best guess for the binary's parameters: that is, we 
determine the best-fit 9\^^ by drawing the normal from 
s to the manifold of approximate waveforms {hAp(6'')}. 
The best-fit parameters 9\^^ satisfy 



{djhAp{9M] 



S-hAp(M =0 



(8) 



Here and below we write 9 for 9^ when the parameter 
index is not required for implicit summations or to iden- 
tify a specific parameter. To first order in the error 
A6'* = 9^^ - 9\^, we can write hAp(6'tr) - hAp(6'bf) « 
-A0^"5jhAp(0bt), so 



S - hAp(6lbf) = n + hGR(6'tr) - hAp(6'tr) + hAp(6'tr) - hAp(6'bf) 
« n + hGR(0tr) - hAp(0tr) - A0^ 9, Hap (^bf ) • 

Defining Tij[9\,f) = ((?.ihAp(6'bt)|9jliAp(^bf)) and inserting Eq. © into Eq. ^ we obtain 

A9' = (t-\9m)Y (^jhAp(0bf)| n) + (r-l(0bf))''' (9jhAp(ebf)| hGR(0tr) - hAp(^tr)) 



(9) 



(10) 
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Thus, to leading order, is the sum of a statistical 
contribution A„0' due to noise and a theoretical contri- 
bution Ath^* due to template inaccuracy, with 



A„0' = (r-i(0M))''(a,hAp(0bf) 



n 



Ath^' = ( r 



d:,hA 



(11) 

hAp(etr) 
(12 



If we happen to know both 6^- and the noise realization 
n, then these two equations can be used to determine 0bf ; 
the experimental situation, however, is that we can deter- 
mine that a certain hAp(0bf) is the best-fit waveform for 
s, but we are left to wonder about the error Ad = 0h[—Sti. 
From this viewpoint, Eq. (fTT|) leads to result ([7]) for the 
distribution of statistical error, which scales as SNR~^; 
by contrast. Equation (jl2p cannot be used as it is written, 
since we do not know ^tr, which for a given 0bf depends 
on the noise. To leading order we can however replace 
hGR(ftr) - hAp(6'tr) by hGR(6'bt) - hAp(6'bf), yielding 



Ath^?^ = ( r 



7bf 



(a,hAP( 



iGRl^bf 



APl^Pbf 



(13) 



which is a noise-independent estimate of theoretical er- 
ror [U, m, [ll]. Moreover, since the signal amplitude 
(and therefore the SNR) appears quadratically in both 
r*-' and in the rightmost inner product of Eq. p3)) . wc 
see that Ath6' is independent of the SNR ^|, so it tends 
to be the limiting factor in extracting source-parameter 
information from the closest, strongest sources. Because 
Ath^* does not decrease with increasing SNR, nor would it 
"average out" if the same event was measured by a large 
number of nearly identical detectors, it can be considered 
a systematic error. This characterization, however, could 
be slightly misleading: while Ath^ is independent of noise 
(to lowest order), it does depend on the true parameter 
values 9tr (e.g., the size and direction errors in the sky 
location depend on the true location). 

Going back to our geometric picture, we see that Eq. 
identifies the signal hGR(0tr) to which hAp(0bf) is 
closest, as it would be given by solving Eq. ^ for ^tr in 
the absence of noise, to leading order in the error Ath^. 
In the next few sections, we develop tools to estimate 
Ath^* that are more accurate than Eq. 



C. A best-fit-by-ODE estimate of theoretical error 

As discussed above, given the observed best-fit wave- 
form hAp(^bf), our goal is to find the true waveform 
h.Qji{9tr) that is best fit by it, in the sense that 



IQR 1(711-^ 



hAp(0 



iGRAC'trJ 



hAp( 



'bf 



(14) 



is minimized, and hence to find the associated parameter 
estimation errors Ath^** = ^bf ~ ^tr- This is a minimiza- 
tion problem, and it could be attacked in many ways (e.g. 



by a straightforward grid search or by simulated anneal- 
ing). Standard methods, however, seem likely to prove 
computationally unaffordable, especially because we wish 
to repeat the minimization for thousands of different 9hf 
in order to explore the distribution of theoretical error 
over parameter space. We have therefore developed an 
alternative minimization scheme that works rather well 
for our particular problem, and that is computationally 
quite efficient. The basic idea is to define a one-parameter 
family of models h\{6) that interpolate smoothly be- 
tween hAp(6') = h.x=o{0) and hGR(6') = hA=i(6'), and 
then to solve the minimization problem on the path from 
A = 0, where we know trivially that 0tr(O) = 6hf, to 
A = 1, where 0tr(l) is the sought-after true-parameter 
vector for the observed 6hf, so that Ath^* = ^bf — ^'tr(l)- 
Namely, each Ou {X) is the solution to 



(hA(0tr(A)) - hAp(0bt) I 9,hAp(0bf)) = 



(15) 



which, once again, implies that (at least locally) hAp(6'bf) 
provides the best fit to hA(6'tr(A)). Taking the derivative 
of Eq. (fTSjl with respect to A yields the ODE 

^(a,hA(0tr(A))|a,hAP(0bf)) = 

- (h^(0tr(A))|9,hAp(^bf)), (16) 

and hence 

^ = -(f-l(0tr(A),0bf))''(hU^?tr(A))|a,hAp(^?M) 



where 



ry (0tr(A),^?bf) = a,hA(0tr(A)) 9jhAp(ebf) 



(17) 



(18) 



Here all partial derivatives are taken with respect to ^bf , 
while h.\{9) = dh\{9)/dX. As anticipated above, the 
initial condition for Eq. P?)) is just 0tr (A)(0) = 6'bf, and 
integrating to A = 1 yields dtr and therefore Ath^ as a 
function of 6'bf. 

Of course, there is an infinite number of ways in which 
we could define the one-parameter family h.x{0) that in- 
terpolates between hAp(6') and hQji{6). Our choice is 
explained in detail in Sec. Ill, but a brief summary can 
be given here. Ignoring modulations due to EISA's mo- 
tion, the phase function '!'(/) has the PN expansion (up 
through 3.5 PN order) 



*(/) = Ax + B + Cx'^/^ 



7 

E 

k=2 



(19) 



where A, B, C and the tpk {k = 2, . . . ,7) are known con- 
stants, and X = {nMfY^^. The signal modulations im- 
posed by EISA's motion depend on the function t(/), 
which represents the instant of time at which the GW 
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frequency sweeps through /, and which has the PN ex- 
pansion 



t{f) ^ D + Ex-^l^ 



k=2 



(20) 



where D, E and the tj, are also known constants. As 
discussed in the introduction, to construct hGR(6') we 
use these full expansions, while for hAp(0) we truncate 
them after fc = 6 (or equivalently, we re-set V'7 and ry 
to 0). For the interpolating family Ha (6'), we make the 
simple replacement 



(21) 



For this choice of \\\(ff), we find that our ODE [Eq. P7)) ] 
leads to a superb fits. Defining the normalized inner 
product (sometimes called match) between two wave- 
forms as 



M(g,k) = 



(g|k) 



(22) 



we generally find (see Sec. IIVI) that the match between 
the true and best-fit approximate waveforms is 



M(hGR(0tr),hAp(f?bf)) > 0.9999. 



(23) 



Incidentally, such high matches mean that it would be 
practically impossible to tell from the detector data alone 
that our theoretical waveforms were inaccurate, since the 
best approximate waveform would be a nearly perfect 
fit to the true waveform, and so would also fit the data 
equally well. 



D. An improved one-step estimate of theoretical 



The ODE method described in the previous section 
suggests an improved version of the linearized estimate 
derived in Sec. Ill Bl We begin by writing each com- 
ponent h'^if) of the interpolating waveform as 



Then 



so we can re-write Eq. pT|) as 



(24) 



(25) 



Eq. {TIJ ate=etr(A) 



djhApiO 



M, 



(26) 

To estimate Ath^*' in a single step (without solving 
the ODE), we approximate A^(6'tr(A)) and *';^(6lti(A)), 



which are functions of the interpolating parameter A, as 
constants: 



^GR 



(^m) 



^AP bf 



) = AA, 



*Uetr(A)) « *GR(^bf) - *Ap(ebf) = A* . 



(27) 



Similarly, while Fy is also a function of A, we find that in 
practice it is approximately constant; thus, to lowest or- 
der in a Taylor-series expansion about A = 0, we replace 
it with 

fy(0tr(A),0bf) ry(0bf) = (a,hAp(0bt)|a,hAp(0bf)) . 

(28) 

With these substitutions the right-hand side of Eq. (|26|) 
becomes constant, and its integration trivial: 



Ath^' - (r-i(0bf))''f [AA + iAA*]e^* 



at e=ebt 



(29) 



Note this is quite similar to the linearized expression ([T 
but with the replacement 



hcR - hAP = Ah ^ [AA + iAA*] e** 



(30) 



Clearly, the two sides of Eq. ([30]) agree up to terms of 
second order in AA and A^f. However, our experience 
shows that Eq. (|29p is more accurate than Eq. (|13p at de- 
termining theoretical errors. Why this is so is explained 
in the next subsection, where we present a second deriva- 
tion of Eq. ([291) . 



E. A second derivation of the improved one-step 
estimate 



As discussed in Sec. IIVI below, we find in practice 
that while Eq. (I13p and our improved one-step estimate 
([29]) agree for very small values of Ath^, Eq. (l29l) re- 
mains a good approximation for a much larger range 
of Ath6'. The basic reason is the following. While in 
practice hAp(0bf) turns out to be extremely close to 
hGR(^tr), with matches exceeding 0.9999, the difference 
hAp(^bf) — hAp(6'tr) (which is the difference of two wave- 
forms from the same family, but evaluated at different 
points) is not well approximated by the first term in its 
Taylor expansion, as assumed in deriving Eq. Q and 
therefore Eqs. P^ : 



hAp(0bf) ~ hAp(0tr) 96 Ae^djhAp{Obf) . 



(31) 



Fortunately, the differences in both the amplitude and 
phase of the waveforms are individually well approxi- 
mated by the linear terms in a Taylor series. 



AAp(0bf) - AAp(0tr) ~ A9^d,AAp{9M) : 
^Ap{Om) - *Ap(etr) ~ A0^dj^Ap{0M) ■ 



(32) 



For replacement PT|) to be reliable, the phase difference 
between the two waveforms should be much less than one 
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radian, A6'^5j *Ap(fi'bf) < 1; however, Eq. is reliable 
as long as A9^A9^didj^Ap{(^hf) ^ Ij which is obviously 
a much less restrictive condition [HHl . 

Given these considerations, we now provide an alter- 
native derivation of our improved one-step formula (j29[) 
along the same lines as our derivation of Eq. (I13|) . We 
begin again with Eq. ([5]), but this time for simplicity we 



neglect the noise term n that determines the statistical 
error. Defining 

(5A = AAp(6'M)-AGR(6ltr), . 
<5* = *Ap(ebf)-*GR(^tr), ^ ^ 

we can rewrite hGR(^tr) — hAp(^bf) as 



AGR(0tr)e'*'=-(^-)-AAp(0bf)e^**-('"' = (AAp(0bf) - <5A)e*(*--(^")-^*) - AAp(0bf)e^**^('"^ (34) 

w -[(5A + i<5*AAp(6ibf)]e'**^(^'^f) . (35) 

In going from Eq. to Eq. (|35p . we have neglected terms of order SAS'if and (5'4'(5'5', which is justified because 
(5A/AAp(^bf) ^ 1 and S'if <SC 1. Next, using Eqs. (j32p and Eqs. (|27p we can rewrite (5A as (5A ~ A0*(?iAAp(^bf 
and similarly J* w A0*aj*Ap(6'bf) + A*. Using d^h = [d,A + iAa,*]e** along with Eq. ([35]), we then have 

hGR(etr) - hAp(ebf) = -A^?*5,hAp(^?bf ) - [AA + zAA*] e'**- , (36) 
and plugging this expression into Eq. ([5]) while setting n 0, we again arrive at 

Ath^* « (r-i(0bf))'' ([AA + *AAvI/]e** a,hAp(0bf)) • (37) 



F. A digression on parameter estimation in the 
face of theoretical uncertainties 



The method for inferring best-fit parameters repre- 
sented by Eq. ^ - simply use our best waveforms, as 
if they were exact, to infer the physical parameters of 
the source - is of course a rather naive way of dealing 
with our theoretical limitations. Since we know that our 
waveforms are not fully accurate, and especially since we 
have some idea of where they fail most badly, we could 
certainly imagine adopting other strategies. For instance, 
since we know that the PN approximation is less reli- 
able for higher values of the PN expansion parameter 
X = (ttM/)^/'^, and therefore at higher frequencies, we 
could somehow modify the inner product ^ to give less 
weight to the higher-frequency part of the waveforms. 
We take a first cut at this approach in Sec. IIVI 

An alternative, Bayesian approach to this problem 
would be to construct a parametrized family of wave- 
forms h.{9,a^) that coincide up through all known PN 
orders, but differ at higher orders. Here the 9^ are again 
the source parameters, while the parametrize our dif- 
ferent models by fixing the unknown higher-PN-order 
terms (so in this context a "model" becomes just a rule to 
derive a waveform h for any given set of physical parame- 
ters 9^.) After assigning an a priori probability p(a^) to 
the different models (refiecting our belief about the de- 
gree to which the model matches GR) , we could marginal- 
ize the Bayesian posterior probability for the 9^ over the 



models, as in 

p{9'\s)^ j p{9'\s-a^)p{a'^)da^ , (38) 

where ^(^'Is; o;-'^) would be obtained by using h.{9'^;a^) 
to build the likelihood. 

Obviously the approach adopted in this paper is much 
less ambitious, but it seems a reasonable first pass at this 
problem; we expect that a more sophisticated approach 
would lead to smaller theoretical errors, so our results 
should provide a rough upper limit to the magnitude of 
Ath^- If the naive error estimate is reassuringly small, 
there may be no need for a more sophisticated treatment. 



III. MODEL OF THE SIGNAL 

To calculate Ath^*' using the formalism developed in 
Secini we need to know the difference Ah = hGR(6'bt) — 
hAp(6'bf) between the exact and approximate waveforms. 
Since calculating Hgr is currently beyond anyone's abil- 
ity, in practice we must content ourselves with estimating 
the rough magnitude of Ath^ from a reasonable estimate 
of Ah. Our basic strategy is to designate the 3.5PN re- 
stricted waveform as hGR, and the corresponding 3PN 
restricted waveform for hAp. We consider two variants 
of this: in one version we write the phasing ^ap through 
3PN order, omitting the 3.5PN terms that go into *gr; 
in the other, we augment ^^'AP with the 3.5PN term of 
lowest order in the symmetric mass ratio rj (which, when 
small, approaches M2/M1). 
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This second version is motivated by the following con- 
siderations. We can clearly improve on the best PN wave- 
forms by employing a hybrid approximation technique 
that expands the Einstein solution as a joint power se- 
ries in both V I c and r\. The lowest-order terms in r\ can 
be obtained to virtually arbitrary order in the PN ex- 
pansion parameter, by means of a perturbative calcula- 
tion for a test particle traveling on a geodesic around a 
BH. This program was begun by Cutler and colleagues 
[2^ and Poisson [2^, for the case of quasi-circular or- 
bits in Schwarzschild, and has since been extended to 
nearl y c ircular, nearly equatorial orbits in Kerr space- 
time |26l [27l . l2q as well as to spinning particles in circu- 
lar, equatorial Kerr orbits [111; for a review, see Sasaki 
and Tagoshi 'sO*]. For the case of circular orbits around 
nonrotating BHs, this PN expansion is now known ana- 
lytically through 5.5PN order [3l|; that is, while binary- 
inspiral waveforms are fully known only through 3.5PN 
order, the lowest-order pieces in 77 are known analytically 
through 5.5PN order. In constructing the best possible 
PN waveforms, it seems appropriate to include the extra 
information derived from BH perturbation theory. We 
shall refer to PN waveforms that have been supplemented 
in this way as "hybrid" waveforms. We do not know 
where this hybrid approach was first clearly spelled out, 
but it is implicit in Kidder, Will, and Wiseman ^S^l, is 
discussed in Damour, Iyer, and Sathyaprakash and 
is used by Canitrot |18j . 

As mentioned in SeclH because this paper represents a 
first pass at the problem of estimating theoretical errors 
for the case of MBHB inspirals, we make some simplifica- 
tions: we work within the restricted PN approximation, 
and we neglect the spin-induced "L x S"' precession of the 
orbital plane, although we do include the lowest-order 
"L • S"' spin-orbit coupling term in the binary-inspiral 
phasing. In addition, we model the LISA response to 
GWs using the "low-frequency approximation" of Cut- 
ler , and the LISA instrument noise using the analytic 
fit of Barack and Cutler [s^j- All these ingredients to 
our computation are briefly summarized in the following 
subsections. 



we use these p° and q"" to define the basis tensors for the 
GW polarization, 

^th = VaVh - qaqb : = PaQb + QaPb ■ (41) 

Within the restricted PN approximation, the gravita- 
tional waveform habit) at the Solar System Barycenter 
can be written as 

habit) = A+it)H+cos<^>it) + A^it)H^^sin<Pit) (42) 

where $(i) is the time-domain PN phasing, and the am- 
plitudes j4+ and Ax are given by 



{A+it),Axit)} = 



|a+,ax| (43) 



rit)D 



with Ml, M2 the two masses, r(i) the instantaneous or- 
bital separation, D the distance to the binary, and 

{a+, ax } = {[! + iL^'na)^] , -2L''na} . (44) 

For the MBHB inspirals considered here, most of the 
LISA SNR accumulates at frequencies / < lOmHz, so 
it is adequate to use the low-frequency approximation to 
the LISA response functions derived by Cutler 0. In 
this approximation, the LISA science data consist of two 
independent (i.e., uncorrelated-noise) channels I and II 
[59| . The LISA channel-I response hiit) is 

hit) = ^A+Fj+((?5,0s>s)cos[$(i) + (^D(i)] 



V3 



AxFl" iOs, c^s,^s) sin [$(t) + ^oit)] (,45) 



where the Doppler phase 



^Dit) = ^Rit)-n 



(46) 



is the difference between the phase of the wavefront at 
the LISA detector [located at i?(t)] and at the Solar Sys- 
tem Barycenter, and where F^^ and F-^ are the detector 
beam-pattern coefficients |35|] 



A. The LISA response to generic 
nonspinning-binary inspiral waveforms 

We describe the binary inspiral as an adiabatic se- 
quence of quasicircular orbits in a fixed plane (since 
we are neglecting spin-induced precession) with nor- 
mal L". Let n"' be the unit vector that points from 
the detector to the source; in terms of the eclip- 
tic colatitude and longitude angles 9s and (t>s, n°- = 
{sin cos ^Sj sin ^5 sin 05, cos ^s} and likewise L° = 
{sin^LCOS<^L,sin^Lsin^L,cos^L}- Now let p° and 
be axes orthogonal to n°, defined by 



e^'^UbPc ; 



(39) 
(40) 



Fj^-''i0s,(l}s,i's) = i(l + cos2 0s) cos 20s cos 2 

=F cos 9s sin 20s sin 2V's • (47) 

Here 9 sit), (f>sit) and ipsit) (without overbars) are the 
source sky-location and polarization angles with respect 
to a frame that rotates with the LISA detector, so they 
vary on a 1-year timescale. Of course, these time- varying 
angles depend on the fixed angles 0s, (j)s, 9l, and 0l, 
although the relations between them are slightly compli- 
cated. We refer the reader to section HI B of Cutler 
for explicit expressions for 9 sit), <j>sit) and ipsit)- 

It is convenient to rewrite the signal ([45l) in the con- 
ventional amplitude-phase form: 



hit) 



y/3MiM2 

rit)D 



Ai it) cos[$ it) -f it) + it)] , (48) 
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where Ai{t) and ipp^i(t) are given by 



fp,i{t) = tan 



a+F+{t) 



(49) 
(50) 



Equations (|45p - ((50|) determine hi{t) completely, while 
the prescription for hii{t) is the same, except that the an- 
tenna patterns F^'^ are replaced by F^"^{9s, (j^S, V's) = 
F+'''{es,(l)s~7r/4,i,s)- 



B. Post- Newtonian templates in the 
stationary-phase approximation 

In this paper we work with the Fourier transform ha{f) 
of the LISA response (where a — 1^ II) , which we approx- 
imate using the stationary phase approximation (SPA) 
i: for / > 0, 



hM) - A4t{f))^M!/'ry' 

X gi(27r(/-/o)to-hA(/)-A(/o)-l-Ao) ^-g-^^ 



■06 



11583231236531 640 tt^ 6848 7 
4694215680 3 21~ 

15335597827 2255 tt^ 1760 6* 12320 5 



3048192 



^1 



76055 o 127825 . 

fl^ n-^ 

1728 ' 

6848 



12 

6848 



9 



1296 



21 



fog(4) : 



21 ' 

/ 77096675 378515 74045 , 
-07 = TT -J^TTTTTT ^ J] 377- r] 



V 254016 1512 



756 



(58) 
(59) 
(60) 



In these expressions, 7 = 0.57721 • • • is the Euler- 
Mascheroni constant, 6 = —1987/3080, and 9 = 
— 11831/9240. Note that ip^ is not listed because it can be 
simply reabsorbed into a re-definition of Aq. In Eq. (j55p 
for 03, /? is the 1.5PN spin-orbit phasing term defined 
in Ref. {(3 is conserved only approximately by the 
1.5PN equations of motion, but in our model we treat 
it as a constant). Because we are not modeling spin- 
orbit precession effects, we have chosen to include only 
the lowest-order spin effects in ^(/); in particular, we 
have not included any spin-spin terms, or any spin-orbit 
terms of order higher than 1.5PN. 
To 3.5PN order, t{f) is given by 



t(/) = io + T(/)-r(/o) 



(61) 



with 



A(/) = *(/) - ^pAKf)) - VdW)) . (52) 



Here Mc is the binary's chirp mass, which is related to 
the total mass M = Mi + M2 and the symmetric mass 
ratio 77 = M1M2/M2 by = Mif/^; also, Aq and Iq 
are the phase and time at the fiducial frequency /o, and 
t = t{f) is the instant at which the GW frequency sweeps 
through the value /, as given below in Eqs. ([6T|l and ([62)l . 
Since ho,{t) is real, ha{-\f\) = K{\f\). 

The main missing ingredient to Eqs. ([?T|) and ([5^ is 
the SPA PN phasing, given up through 3.5PN order in 
Refs. Hill: 



vI/(/) = - (SvrMj) [1 + ^ i^kx'' + ^ i^kx'' log(x) 



k=2 



(53) 

with the PN expansion parameter x{f) = (ttA//)^/'^, and 



20 /743 11 
^ Tl 336 + T'^ 



(54) 
(55) 



i/,3 = -i67r-h4/3, 

/3058673 5429 617 2 , 



05 



V 1016064 1008 ' 144 
38645 65 



252 3 ' ' 



(57) 



r{f) = 



S^/)-«/3A4-5/3 

r 7 6 

X l + ^Tkx'' + ^fkX^\0g{x) 



k=2 



k=5 



(62) 



where Tk = ^Vfe for /c < 5, T5 = -|05, tq = -■|(V'6 + 
■06), 7=6 = -1-06, and T7 = -|V'7- 

We truncate the waveform of Eq. (|5ip above a fre- 
quency corresponding to either r = 6M or (to explore 
the theoretical error due to the last phase of inspiral) 
r = 9M, and below an initial frequency chosen by solving 
Eq. ([S^ to enforce a signal of a certain fiducial length 
(typically one year). Thus, our signal model includes 
just the inspiral waveform, and not the final merger and 
ringdown. This is in part because of our current igno- 
rance about the final merger: NR computations are now 
producing robust merger waveforms for nonrotating BHs 
|36|, and the case of generic spins is beginning to be ex- 
plored [13, m, lio, Hlf, but the relevant parameter 
space is very large, and much work is still needed. In 
addition, in the particular problem considered in this pa- 
per, most of the interesting information about binaries 
(and especially about their sky locations) is presumably 
encoded in their phase evolution during the slow inspi- 
ral, so one may hope that neglecting the information in 
the merger and ringdown does not strongly impact pa- 
rameter estimation, even in cases where the latter phases 
dominate the total SNR. Establishing whether this is true 
certainly deserves careful investigation in the future. 
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Overall, ha{f) depends on ten physical parameters: 
the mass parameters M^. and rj, the spin-orbit param- 
eter (3, the phase Aq and time to at the fiducial fre- 
quency /o , the ecliptic sky-position angles 9s and (jjs and 
orbital-angular- momentum angles 0l and (j)L (which to- 
gether determine the LIS A- frame angles 9s, (f>s, and ips), 
and the distance D. Note that for simplicity we treat 
the background spacetime as flat, rather than Friedman- 
Lemaitre-Robertson-Walker. Accounting for cosmologi- 
cal effects simply requires the translation Mi Mi{l+z) 
and D ^ Dl, where is the luminosity distance [i^ . 



The fiducial, approximated, and interpolating 
waveform families 



Our fiducial (hcR), approximated (Hap), and inter- 
polating (Iia) template families are based on the LISA 
response model and SPA PN waveforms outlined in Sees. 
ImXl and lllTB] 

• To build liGR, we evaluate and t{f ) from Eqs. 
([55)) - (p^ . including all terms given there (through 
3.5PN). 

• We build two distinct versions of Hap: a straight 
version that is the same as hcR, except that we set 
^/^r and ry to zero by hand; and a hybrid version 
where we set V'7 and t-j to their 0{Tf) terms alone. 



^ /77096675\ 
^ ^ V 254016 J ' 



(63) 



and Ty —(2/5) x (the new ip-j). 

(Recall the motivation for the hybrid version: we 
assume that terms of ipj and Tj of lowest order in 
T] are known from BH perturbation theory, so hAP 
must converge to hcR as 77 ^ 0.) 

• We also use two versions of the interpolating family 
h^. The straight version is the same as Hgr, except 
that we modulate the strength of the 3.5PN terms 
with the replacement ^7 X'ipj and T7 — ^ Ary; 
clearly this family coincides with the straight Hap 
when A = 0, and with hcR when A = 1. 

The hybrid version is obtained with the replace- 
ment 



■07 



/77096675\ 


+ A 






H 


V 254016 J 







/ 378515 
1512 



74045 
756 



(64) 

and T7 —(2/5) x (the new V'y), so that again 
h> = hAp for A = and = hep, for A = 1. 



D. LISA noise model 

We assume that the LISA noises in channels I and 
II are uncorrelated, so the inner product between the 



signals g = {gi{t) , gu{t)) and k = {ki{t) , ku{t)) is just 



(g|k) =4Re J2 / 

0=1,11-^0 



■df. (65) 



We also also assume that S\{f) and S]^{f) are exactly 
equal, and given by the fitting function Sh{f) of Barack 
and Cutler (Note however that the Sh{f) used 

here is actually 3/20 times the expression quoted in Ref. 
[3^ . This is because Ref. 0| uses the sky-averaging 
convention, but this paper does not.) The LISA noise 
has three main components: instrument noise, confusion 
noise from short-period galactic white-dwarf (WD) bina- 
ries, and confusion noise from extragalactic WD binaries. 
(We neglect confusion noise from unresolved extreme- 
mass-ratio inspirals, since its magnitude is quite uncer- 
tain, and since it is likely to be dominated by either in- 
strument noise or WD confusion noise at all frequencies 

For the LISA instrument noise we adopt the fitting 
function of Finn and Thorne [43[ , which is based on the 
noise budgets specified in the LISA Pre-Phase A Report 

ST\f) - (9.18 X IQ-^^if/Uz)-^ + 1.59 x lO'-^^ 

+ 9.18 X 10"3^(//Hz)2)hz"\ (66) 



Next we turn to confusion noise from galactic and extra- 
galactic WD binaries (GWDs and EWDs, respectively). 
Any isotropic background of individually unresolvable 
GW sources represents (for the purpose of anal yzin g 
other sources) a noise source with spectral density [45| 



S 



57r^ d(log/) 



(67) 



estimates of dpcw/rflog/ from the GWD and EWD 
backgrounds [4a, |43| yield the spectral densities 



cGWD 



CEWD 



(/) = 2.1 X 10-^^ ( I- 
(/) = 4.2 X 10-4^ ( ^ 



-7/3 



Hz" 



-7/3 



Hz' 



(68) 
(69) 



Essentially none of the EWDs can be individually re- 
solved and fit out, so to a first approximation one can 
treat Sj^'^^lf) as if it were just another source of in- 
strument noise. The situation is different for the GWDs. 
While this background is actually larger than the LISA 
instrument noise between ~ 10~* and 10~^ Hz, at fre- 
quencies / > 2 X 10"'^ Hz the galactic sources are suffi- 
ciently sparse, in frequency space, that one expects to be 
able to fit them out of the data. This process effectively 
reduces the amount of information about any other sig- 
nal h in the data — this is because the piece of h that lies 
in the tangent space to the signal manifold that describes 
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the resolvable WD binaries gets fitted out as part of the 
WD background [48| . 

Consider a frequency band of width A/, much smaller 
than the LISA bandwidth, but much larger than l/Ttot, 
where Ttot is LISA's data-collecting lifetime. The num- 
ber of independent data points collected in that band 
is 4TtotA/, where the factor 4 is the product of 2 data 
channels (I and II) times 2 real numbers (or one complex 
number) per discrete frequency bin. The dimensionality 
of the subspace that is given up to the WDB background 
is 7{dN/df)Af, where dN/df is the density of GWDs in 
frequency space and 7 is the number of parameters re- 
quired to specify each GWD. Thus the fraction of the in- 
formation in h that is effectively lost to the background 
is (J / A) dN / df T^^l , and ShJJ) is effectively multiplied 



by [l-(7/4)diV/d/rt; 



This simplistic counting 



clearly breaks down as {7 / A) dN / df T^l 1, since the 
total effective noise should not exceed the sum of all the 
individual noises. Lacking a more sophisticated treat- 
ment of this transition, we adopt the following fitting 
function as our (admittedly crude) estimate of LISA's 
total effective noise density Sf^{f): 

Sfif) = min{ [Sr\f) + 5fw°(/)] exp{KT~'dN/df), 
For dN/df we adopt the estimate 



(70) 



dN 



= 2 X IQ-'^Hz" 



IHz 



11/3 



(71) 



In the paragraph above Eq. (I70p we have given an ar- 
gument suggesting that n ~ 7 /A (assuming that both 
channels I and II are operational), but this really rep- 
resents the best one could do. Partly to be consistent 
with earlier papers, we account for suboptimal WD sub- 
traction by adopting a more conservative estimate for k 
given by nT^l = 1.5/yr. Because dN/df falls so steeply 
with increasing /, had we (very optimistically) taken 
uTt'l = 7/(20yr), the steep drop off in Sf{f) at / - 2 
mHz would merely have been shifted a factor ^ 1.5 to- 
wards lower /. 



waveforms and LISA response change very smoothly with 
these parameters, so this approximation is adequate). In 
addition, the following details about our implementation 
are worth pointing out: 

• The PN-interpolation parameter A enters the LISA- 
response expressions of Sec. lIII Al through t{f) [Eqs. 
(|6ip and (|62p ]: we include this dependence when 
computing numerical derivatives. 



In principle, the derivative dh/dlog Mc includes 
terms of the form {dh/dFiji) x (dFiji/dt) x 
{dt/d\ogMc). However, these contributions are 
quite time-consuming to compute, and they are also 
^ 10'^ times smaller than the dominant terms of 
the form {dh/d'i') x (dip/dlogMc)- Therefore we 
neglect the former, as well as similar terms in the 
derivatives dh/dr] and dh/d(3. 



• We choose the Nyquist frequency and frequency 
spacing used to represent waveforms in the fre- 
quency domain in such a way to avoid both fre- 
quency aliasing and time wrapping. 



• The artificial truncation of waveforms at GW fre- 
quencies corresponding to r = 6M (or 9M) induces 
a mildly pathological behavior in the inner products 
of signals from systems with different total masses; 
we deal with this problem by truncating all wave- 
forms at the r = 6M (or 9M) GW frequency of 
the hAp(^bf) for which we are estimating theoret- 
ical error. This approximation is justified because 
the frequency-domain truncation does not repre- 
sent any physical feature of the true waveforms, 
but stems instead from our ignorance about the 
late stages of their inspiral. Furthermore, the trun- 
cation does not affect the information carried by 
the signals at lower frequencies, and while it does 
exclude some high-frequency information, it is in- 
formation that we are ill-equipped to interpret any- 
way. 



IV. RESULTS 

To test our formalism numerically, we have developed 
MATLAB code to do the following: a) generate the gravi- 
tational waveforms, LISA response, and LISA noise spec- 
trum described in Sec. IIIIl b) compute noise inner prod- 
ucts according to Eq. ([3]); and c) implement the ODE 
[Eqs. ((ni)-(ITHl)] and improved one-step [Eq. formu- 
las for the theoretical error. We use analytical derivatives 
of the signals with respect to all parameters except the 
PN-interpolating parameter A and the sky-position and 
L angles, for which we use central finite differences with 
displacements of 10~^ and 2tt x 10~^, respectively (the 



• Last, we set the fiducial frequency fo to 90% of 
the truncation frequency; since most of the total 
SNR for the systems studied in this paper comes 
from the last stages of inspiral, this has the effect 
of reducing the covariance in the joint estimation of 
Aq and to [see Eq. ([51]) ] , and therefore of increasing 
the robustness of our framework (and especially of 
the one-step formula). 

In the remainder of this section, we test the accuracy 
of the improved one-step formula against the solution 
of the ODE, and we survey theoretical errors for MBHB 
systems with component masses between lO** and IO^Mq. 
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A. Comparison of the ODE and improved one-step 
formulas 

To evaluate the accuracy of the improved one-step for- 
mula against the ODE formula, we select a representa- 
tive system with Mi = 0.5 x IO^Mq, A'h = IO^Mq, 
and /3 = 0, and we compute expected theoretical errors 
for 600 random choices of the sky-position and L angles, 
uniformly distributed on their respective spheres. (Recall 
also that the theoretical errors are independent of SNR, 
so there is no need to specify the latter.) All signals are 
created with a duration of one year (i.e., they begin a year 
before the system reaches the truncation frequency) ; the 
LISA orbits are chosen according the conventions of Ref. 
0, and time-shifted so that the LISA position and ori- 
entation are the same for all binaries at the time when 
they achieve their respective fiducial frequency /q. 

The ODE and one-step errors are compared in the scat- 
ter plots of Fig.m where the clusters labeled "straight," 
"trunc," and "hybrid" refer, respectively, to errors com- 
puted for our original waveforms truncated at r = 6A/, 
for waveforms truncated at r = 9Af, and for waveforms 
with the hybrid phasing described in Sec. lIII Cl fand again 
truncated at r = 6M) . The margins of each scatter plot 
display the single-variable distribution of the two errors 
over the populations of sky positions and L orientations. 
We see that for the "intrinsic" parameters Mc, 77, /?, for 
the sky-location angles, and for the time t^, the one- 
step formula provides a reasonable approximation for the 
ODE result, which improves in going from the straight 
signal model to the truncated model, and (even more so) 
to the hybrid model. 

In fact, the one-step errors appear to deviate from 
the ODE results by almost-constant multiplicative fac- 
tors (approximated in the figure by least-squares-fitting 
the clusters to lines), which indicates that the error in- 
curred in reducing the ODE to a single step is roughly 
independent of sky-position and L angles. Unfortunately, 
at present we do not have a way to predict the corrective 
factors from the parameters of the system. The behavior 
of the one-step formula is more erratic for the remaining 
parameters (the L angles 6l and (pLi the phase Ag, and 
the amplitude A), with many outliers that are several 
times the maximum reported ODE error, which makes 
it hard to prepare informative scatter plots. However, 
the median errors (see the table within Fig. [J) are again 
close, and the outliers seem to be due to numerical error 
in the computation of F^^ for quasi-singular F's. 

As a confirmation of the accuracy of the ODE method, 
we found the "final" match Af(hAp(^?bf )7 hGR(6'tr)) [see 
Eq. ^] to be > 0.9999 for all sets of angles, and for 
all three waveform variants; final matches were consid- 
erably lower with the one-step formula (except for the 
hybrid waveforms, where they were > 0.99 for 93% of 
the Monte Carlo population); and the "initial" matches 
M(hAp(0bf), hGR(6'bt)) were always lower than 0.50. 

We conclude that the improved one-step formula yields 
at least the correct order of magnitude for theoretical er- 



ror, and usually does much better than that. In partic- 
ular, it approximates well the theoretical errors for the 
intrinsic parameters and for the sky location, which is 
of special interest in view of the possibility of searching 
for electromagnetic counterparts to LISA MBHB sources. 
The one-step formula is attractive both because it is easy 
to implement and because it runs much faster than the 
ODE, by a factor equal approximately to the typical 
number of integration steps (in our work we found that 
^ 40 steps were usually sufficient to achieve several digits 
of accuracy in 0tr)- It is therefore appropriate to use the 
improved one-step formula in rapid Monte Carlo surveys 
over large parameter ranges, such as the survey discussed 
in the next section. 



B. Monte Carlo survey of theoretical errors 

We survey theoretical errors for eight representative 
MBHB systems defined by the component mass combina- 
tions {M1+M2) = {10^ + 10^)Mq, (10'*-hl06)MQ, (10^-1- 
106)Mo, (10'^ + 10^)Mq, (10*-hlO^)M^(105-HlO^)Mo, 
(10^ + 1O^)M0, and (10^ + 10^)Mo [g^ and by /3 = 0. 
To do so, we compute improved-one-step theoretical er- 
rors for the straight, truncated, and hybrid signal mod- 
els, for 600 random choices of sky-position and L angles 
uniformly distributed on the sky. In addition, we com- 
pute expected statistical errors from the Fisher-matrix 
formula [Eq. ([7])], using 3.5PN waveforms truncated at 
r = 6M, and setting the SNR = 1000. SNRs of this 
magnitude are indeed expected for the best LISA MBHB 
detections '3] , so they provide a reasonable benchmark 
to gauge the magnitude of theoretical errors. 

To visualize the relation between theoretical and statis- 
tical errors, consider Fig. [T] the thick dot at A log Mc ~ 
A77 — marks the location ^bf of the assumed best-fit pa- 
rameters log Mc and 77, while the small solid and empty 
dots show the location 9tr of the true system parame- 
ters (computed for this figure with the ODE formula and 
the hybrid signal model), for two representative sky loca- 
tions and L orientations, described at the bottom left of 
the figure. The ellipses around the dots enclose Icr prob- 
ability regions for 0tr in the presence of noise-induced 
statistical error A0„ (computed from the Fisher matrix, 
and approximated as independent from A^th as in Sec. 
Ill Bp . Recall that for a given 6'bt, A^th is deterministic 
(although it is unknown in practice, because we do not 
have access to the true Hgr), while A0„ is a random 
variable with distribution determined by the properties 
of the waveforms and of instrument noise. 

In Tab. [T] we show the median of the absolute theoret- 
ical errors over the populations of 0s, <t>S, ^Lj and 4>l, 
for each representative mass combination, and for our 
three signal models. The first row in each group gives 
the median la statistical errors. In the last column, A is 
defined as the combination Mc^^ /D. [Median absolute 
errors are more robust indicators of typical errors than 
the error variances, which can be corrupted (especially 
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FIG. 1: Theoretical errors (dots) and la statistical errors (el- 
lipses) for several MBHB mass configurations, and for two 
representative sets of sky-position and L angles (near the 
ecliptic plane — > solid dots and ellipses; closer to the pole 
—> empty dots and dashed ellipses). For all systems /3 = 
and SNR = 1000. Strictly speaking, this plot displays minus 
the theoretical error, since the dots mark the location of 9ti 
relative to ^bf (the thick dot at the origin). Statistical errors 
become smaller for more asymmetric mass configurations, and 
for smaller total masses. The most asymmetric mass config- 
uration, (10* -I- 1O^)M0, displays the largest spread of 6ti. 



for A9l, A(^l, AAq, and A log A) by the few outliers 
produced by the improved-one-step formula. The single- 
parameter la statistical errors correspond in Fig. [T] to 
half the projection of the ellipses on the horizontal and 
vertical axes.] 

Several trends are apparent from Tab. IH First, theo- 
retical errors decrease (as expected) when we move from 
the straight signal models (truncated at r = 6M) to the 
early-truncation (at r = 9M) models. The median im- 
provement is only a factor ^ 1.5-2.5 for the intrinsic 
parameters Mc, rj, and (3, but ~ 2-9 for the angles 0s, 
Sl, and (pL- Theoretical errors are even smaller for 
our hybrid signal model, where Hap tends to hcR in the 
limit r] 0. The improvement is most dramatic (factors 
of 10^-10^) for mass combinations that actually have very 
smah mass ratios, such as {M1 + M2) = (1O^-|-1O'5)M0 or 
(10^ -I- 10'^)Mq. But the improvement over the straight 
model is substantial (factors of ~ 10) even for the equal- 
mass cases. 

While our estimated theoretical errors are generally 
small on an absolute scale, they are of comparable mag- 
nitude or larger than the SNR — 1000 statistical er- 
rors, especially for the intrinsic parameters Mc, rj, and 
(3, and for Aq and io [with one exception for the hybrid 
model at the highest mass ratio, (10* -I- 1O'')M0]. For 
the sky-position angles Os and 0s 1 theoretical and sta- 
tistical errors are comparable when we use the straight 
and early-truncation signal models; however, Ath^s and 



Ath0s are always comfortably smaller than the statisti- 
cal errors when we use the hybrid model. The results for 
9l, <j)L, and log^ are quite similar to those for ^5 and 
05 (although, as mentioned above, our one-step formula 
is less reliable for those parameters). 



SUMMARY AND FUTURE WORK 



The results of Sec. IIVBI suggest that for hybrid PN 
waveforms, theoretical errors due to waveform inaccu- 
racy will not be large enough to significantly degrade the 
available science from MBHBs. Nevertheless theoretical 
errors could well be the limiting factor in determining 
source parameters for the strongest sources. This is espe- 
cially true for the intrinsic parameters such as the masses 
and spins; theoretical errors appear relatively more be- 
nign (compared to statistical errors) for the sky posi- 
tion angles. [Also, to enable searches for electromagnetic 
counterparts to the inspiral GW signals, the sky position 
will have to be determined several days before the end of 
the inspiral, to allow for the intermittent transmission of 
LISA data to Earth, and to provide time to set up elec- 
tromagnetic observations I50|- This will reduce the SNR 
available for detection (and therefore increase statistical 
error) , but it will also reduce the theoretical error due to 
unknown high-PN-order terms, which are suppressed by 
fractional powers of the orbital frequency.] 

The indications of this paper are tentative, in that 
the waveform models adopted in this paper are ex- 
tremely simplified; several additional known effects, in- 
cluding amplitude modulations due to dynamically evolv- 
ing spin-orbit and spin-spin couplings (see, e.g., Ref. 
51 1), significant orbital eccentricity @, higher-PN-order 
GW harmonics ^5^, and post-inspiral parts of the wave- 
form (i.e., merger and ringdown), modify the signal-space 
geometries of the template and true waveform manifolds 
in ways that are beneficial and ways that are damaging 
for parameter estimation; the net effect needs to be com- 
puted. Improved surveys of theoretical error that include 
the waveform features listed above can be performed with 
the mathematical tools described in Sec.|TT]of this paper. 
This will be very important in establishing whether fur- 
ther work in PN and NR calculations is needed for LISA 
to reach its full scientific potential (ultimately limited by 
the source-confusion and instrument noise profiles, and 
not by the development of GR theory). 

The fact that the match Af (hAp(6'bf), hGR(6'tr)) turns 
out to be > 0.9999 for all three approximate waveforms 
is a rather striking result. If this continues to hold for 
more realistic waveform models, it will mean that it may 
be quite hard (for typical SNRs) to judge, from goodness 
of fit to the data, how accurate one's theoretical template 
waveforms really are. 

There are other applications of those methods. For 
instance, almost the same problem studied for LISA in 
this paper arises in the context of ground-based GW as- 
tronomy, in estimating parameters from the inspirals and 
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TABLE I: Median absolute theoretical errors in our survey of 600 sets of sky-position and L angles for eight representative 
mass combinations (with /3 = 0). The first line in each group gives the Icr single-parameter statistical error at SNR = 1000; the 
second to fourth lines give the median absolute theoretical errors for the straight, early-truncation, and hybrid signal models. 





A log Mc 


Arj/rj 


A/3 


AOs 


A0S 




A0i 


AAo 


foAto 


A log A 


nn^ + 10^~lMr^ Stat 
sys. 

(trunc.) sys. 
(hybrid) sys. 


9.98 10"'' 
4.4710"* 
2.8410"* 
5.86 10"^ 


2.29 10"^ 
1.89 10"^ 
1.40 10"^ 
2.46 10"^ 


1.76 10"^ 
1.18 

8.30 10"^ 
1.53 10"* 


1.96 10"^ 
5.46 10"^ 
2.51 10"^ 
7.30 10"^ 


2.44 10"^ 
8.9410"^ 
3.9710"^ 
1.22 10"^ 


3.04 10"^ 
1.0710"* 
4.6410"^ 
1.45 10"^ 


3.96 10"^ 
1.55 10"* 
6.49 10"^ 
2.11 10"^ 


9.55 10"^ 
3.17 

7.84 10"* 
4.2710"^ 


1.23 10"* 
1.58 

9.73 10"^ 
2.02 10"^ 


2.45 10"^ 
9.41 10"^ 
3.8810"^ 
1.2610"^ 


no'' + ^Vl^^Mr^ stat 

sys. 

(trunc.) sys. 
(hybrid) sys. 


1.38 10"^ 
3.59 10"^ 
1.68 10"^ 
2.00 10"^ 


5.82 10"* 
3.20 10"^ 
2.02 10"^ 
1.7710"^ 


5.07 10"^ 

2.04 

1.13 

1.13 10"^ 


1.49 10"^ 
1.39 10"^ 
1.62 10"^ 
8.21 10"* 


1.74 10"^ 
1.80 10"^ 
2.09 10"^ 
1.08 10"^ 


1.87 10"^ 

1.88 10"* 
2.42 10"^ 
1.12 10"^ 


2.66 10"^ 
2.76 10"^ 
3.32 10"^ 
1.60 10"^ 


1.57 10"* 
1.52 

1.52 10"* 
8.99 10"^ 


2.07 10"^ 
2.89 10* 
1.79 10* 
1.59 10"* 


1.67 10"^ 
1.5310"* 
2.23 10"^ 
8.88 10"* 


(10^ + 1O'')M0 stat. 

sys. 

(trunc.) sys. 
(hybrid) sys. 


3.27 10"^ 
1.90 10"^ 
8.36 10"* 
8.70 10"^ 


2.53 10"^ 
2.83 10"^ 
1.68 10"^ 
1.29 10"^ 


2.04 10"^ 
1.74 

9.05 10"^ 
7.92 10"^ 


2.53 10"^ 
6.01 10"^ 
8.42 10"^ 
2.82 10"^ 


3.05 10"^ 
7.21 10"^ 
1.01 10"^ 
3.46 10"^ 


3.66 10"^ 
1.16 10"^ 
1.29 10"^ 
5.40 10"^ 


4.65 10"^ 
1.23 10"^ 
1.68 10"^ 
5.85 10"^ 


2.64 10"^ 
6.54 10"* 
6.2710"^ 
3.1310"^ 


2.81 10"^ 

3.80 

2.40 

1.71 10"^ 


2.89 10"^ 
8.33 10"^ 
1.01 10"^ 
3.9610"^ 
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sys. 

(trunc.) sys. 
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3.32 10"^ 
1.89 10"^ 
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1.21 

2.6410"^ 
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2.15 10"^ 
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2.31 10"^ 
1.89 10"^ 
2.86 10"^ 
2.26 10"^ 


2.89 10"^ 
3.69 10"^ 
4.89 10"^ 
4.50 10"^ 


3.77 10"^ 
4.39 10"^ 
5.86 10"^ 
5.26 10"^ 


1.22 10"^ 
1.49 10"* 
1.38 10"^ 
1.76 10"^ 


1.28 10"^ 
1.46 

9.36 10"* 
1.64 10"^ 


2.50 10"^ 
3.2810"^ 
3.93 10"^ 
3.9410"^ 


(10'' + 1O'^)M0 stat. 

sys. 

(trunc.) sys. 
(hybrid) sys. 


3.56 10"^ 
1.4410"^ 
7.22 10"^ 
6.65 10"*^ 


6.66 10"* 
5.51 10"i 
3.82 10"^ 
2.55 10"* 


6.80 10"^ 

4.33 

2.70 

2.00 10"^ 


1.46 10"^ 
1.46 10"^ 
5.2410"^ 
9.53 10"® 


1.75 10"^ 
1.63 10"^ 
8.2710"^ 
1.12 10"^ 


2.40 10"^ 
2.63 10"^ 
1.03 10"^ 
1.30 lO"'' 


3.37 10"^ 
3.82 10"^ 
1.64 10"^ 
1.68 10"^ 


4.81 10"^ 
8.96 10"2 
3.46 10"^ 
3.88 10"^ 


4.68 10"^ 
2.62 10^ 
1.54 10^ 
1.21 10"* 


2.47 10"'^ 
1.81 10"^ 
4.64 10"^ 
9.43 10"® 


(10^ + 1O^)M0 stat. 

sys. 

(trunc.) sys. 
(hybrid) sys. 


3.51 10"^ 
6.5410"^ 
3.36 10"^ 
3.3710"^ 


1.04 10"^ 
3.9710"^ 
2.7410"^ 
2.0410"^ 


9.60 10"^ 

2.78 

1.73 

1.43 10"^ 


2.29 10"^ 
9.62 10"^ 
2.39 10"^ 
6.49 10"^ 


2.77 10"^ 
8.22 10"^ 
2.4410"^ 
5.12 10"^ 


3.93 10"^ 
1.53 10"^ 
3.66 10"^ 
7.8710"^ 


5.08 10"^ 
2.08 10"^ 
5.01 10"^ 
1.12 10"* 


6.89 10"^ 
3.39 10"^ 
9.48 10"^ 
1.85 10"* 


1.40 10"^ 
2.83 10* 
1.70 10* 
1.45 10"* 


3.57 10"^ 
9.23 10"^ 
2.4310"^ 
6.02 10"^ 


(10® + 1O'^)M0 stat. 

sys. 

(trunc.) sys. 
(hybrid) sys. 


6.66 10"^ 
3.45 10"^ 
1.7410"^ 
1.49 10"* 


3.23 10"^ 
3.41 10"* 
2.28 10"* 
1.46 10"^ 


2.82 10"^ 

2.29 

1.38 

9.85 10"^ 


4.69 10"^ 
6.38 10"^ 
1.4410"^ 
3.08 10"* 


5.5410"^ 
5.06 10"^ 
1.12 10"^ 
2.42 10"* 


6.96 10"^ 
9.46 10"^ 
2.29 10"^ 
4.85 10"* 


9.79 10"^ 
1.32 10"^ 
2.58 10"^ 
5.7710"* 


1.11 10"^ 
1.50 10"^ 
3.46 10"^ 
7.63 10"* 


6.9710"^ 

3.75 

2.29 

1.59 10"^ 


6.46 10"^ 
7.25 10"^ 
1.6310"^ 
3.56 10"* 


(10^ -f 1O^)M0 stat. 

sys. 

(trunc.) sys. 
(hybrid) sys. 


1.65 10"* 
3.13 10"^ 
1.50 10"^ 
3.29 10"* 


1.08 10"^ 
4.35 10"^ 
2.8710"* 
4.53 10"^ 


1.01 10"* 

3.31 

2.02 

3.45 10"* 


6.8110"^ 
4.9410"^ 
1.05 10"^ 
5.8110"* 


8.23 10"^ 
3.28 10"^ 
7.33 10"* 
3.8710"* 


1.18 10"^ 
7.10 10"^ 
1.83 10"^ 
7.9710"* 


1.59 10"^ 
1.0710"^ 
2.2710"^ 
1.2710"^ 


1.5710"^ 
1.05 10"^ 
2.4710"^ 
1.34 10"^ 


5.53 10"^ 
1.43 

8.65 10"* 
1.4710"* 


9.8910"^ 
6.1910"^ 
1.3910"^ 
7.4710"* 



mergers of binaries of neutron stars and/or stellar-mass 
BHs. Also, there are many cases in GW observations 
where theoretical approximations are made to simplify 
calculations (again effectively modifying the "search" 
template family with the respect to the true measured 
signals), and the tools of this paper could be used to 
investigate the impact of these approximations on the 
extraction of source information from observations. One 
such approximation, recently critiqued by Grishchuk [s^ , 
is the long- wavelength formula for the LIGO and VIRGO 
GW response functions, which introduces distortions at 
frequencies of the order of the inverse light-travel time 
across the interferometer arms. Another application is 
to LISA searches for extreme-mass-ratio inspirals (EM- 
RIs), where 77 is typically ~ 10"^. A great deal of effort 



has been expended, and impressive progress made, in cal- 
culating EMRI waveforms to first-order in an expansion 
in r/ [54]. How important is it for theorists to calculate 
EMRI waveforms through order 77^ [1^1? We plan to ad- 
dress that question in a future paper. 
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FIG. 2: Scatter plots of the ODE theoretical errors (horizontal axes) against the improved-one-step theoretical errors, for MBHB waveforms with Mi = 0.5 x 10® Mq, 
M2 = 10'' Mq, and /3 = 0, over populations of 600 sets of sky-position and L angles. The clusters labeled "straight," "trunc," and "hybrid" refer to errors for the 
three waveform variants considered in this paper (truncated at r = 6M, at r = 9M, and with hybrid phasing). The margins of each plot show the one- variable 
distribution of the two types of error, with vertical bars marking the location of the median, and of the 5% and 95% quantiles (except where stated differently); these 
plots were normalized arbitrarily, but consistently for each parameter. The lines drawn through the "straight" and "trunc." clusters represent least-squares fits to 
A^ono-stcp(ASoDE), with the correspondent linear coefflcient indicated nearby (by comparison, the dashed lines are the locus of A6'onc-stop = A^ode). Scatter plots are 
not shown for parameters 9l, <f)L, Ao, and A for which the one-step formula can be farther off the ODE prediction (although medians are still close, as shown in the 
in-figure table). 



